1 Data import

The data is included within this reproducible report and can be downloaded here.

1.1 Data wrangling

# data frames with sum scales
data_scales <- data_rbt%>%
  mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
         Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3, 
                                   exp_4, exp_5, exp_6), na.rm = T),
         Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
         Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
         TSM = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4), na.rm = T))

data_scales_lables <- data_rbt%>%
  mutate(Treatment = factor(treat, levels = c("GB", "CC", "CB")),
         Treatment = fct_recode(Treatment,
                                `Colored badges` = "CB",
                                `Control Condition` = "CC",
                                `Greyed out badges` = "GB"),
         Experitse = rowMeans(data.frame(exp_1, exp_2, exp_3, 
                                   exp_4, exp_5, exp_6), na.rm = T),
         Integrity = rowMeans(data.frame(int_1, int_2, int_3, int_4), na.rm = T),
         Benevolence = rowMeans(data.frame(ben_1, ben_2, ben_3, ben_4), na.rm = T),
        `Topic specific multiplism` = rowMeans(data.frame(tsm_1, tsm_2, tsm_3, tsm_4), 
                                               na.rm = T))

data_scales_wide <- data_scales%>%
  select(Experitse, Integrity, Benevolence,
         TSM, Treatment, session)%>%
  gather(Variable, Value, -session, -Treatment)%>%
  mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
  select(-Variable, -Treatment)%>%
  spread(Variable2, Value)

data_scales_lables_wide <- data_scales_lables%>%
  select(Experitse, Integrity, Benevolence,
         `Topic specific multiplism`, Treatment, session)%>%
  gather(Variable, Value, -session, -Treatment)%>%
  mutate(Variable2 = paste(Variable, Treatment, sep = "_"))%>%
  select(-Variable, -Treatment)%>%
  spread(Variable2, Value)

# Data of treatment check
data_tch <- data_rbt %>% 
  select(starts_with("tch"), meas_rep, treat)

data_tch_n <- data_tch%>%
  mutate_all(function(x) ifelse(x == -999, NA, x)) %>% 
  filter(rowSums(is.na(data.frame(tch_1, tch_2, tch_3, tch_4, tch_5))) != 5)

PID_list <- data_rbt$session %>% unique()

2 Sample

2.1 Sample size

length(PID_list)
## [1] 270

2.2 Skipped after first measurement

sum(is.na(data_rbt%>%
            filter(meas_rep == 2)%>%
            pull(exp_1)))
## [1] 0

2.3 Demographics

data_rbt %>% 
  select(age, semester, sex) %>% 
  skim(.)
Data summary
Name Piped data
Number of rows 527
Number of columns 3
_______________________
Column type frequency:
numeric 3
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
age 21 0.96 22.89 2.95 18 20 22 25 33 ▇▆▅▁▁
semester 21 0.96 5.86 3.68 1 3 5 9 19 ▇▅▅▁▁
sex 21 0.96 1.71 0.47 1 1 2 2 3 ▃▁▇▁▁

2.3.1 Count on sex

data_rbt$sex %>% 
  table(.)
## .
##   1   2   3 
## 150 352   4

3 Psychometric properties of the measurements

3.1 Descriptive overview

skewn <- function(x) DescTools::Skew(x, na.rm = T)
kurto <- function(x) DescTools::Kurt(x, na.rm = T)
maxabszvalue <- function(x) max(abs(scale(na.omit(x))))

my_skim <- skim_with(numeric = sfl(skewn, kurto, maxabszvalue))

data_scales%>%
  my_skim(.)
Data summary
Name Piped data
Number of rows 527
Number of columns 34
_______________________
Column type frequency:
character 2
factor 1
numeric 31
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
session 0 1 64 64 0 270 0
treat 0 1 2 2 0 3 0

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
Treatment 0 1 FALSE 3 GB: 184, CC: 173, CB: 170

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist skewn kurto maxabszvalue
exp_1 0 1.00 5.43 1.29 1 5.00 6.00 6.0 7 ▁▁▂▃▇ -0.87 0.40 3.42
exp_2 0 1.00 5.44 1.15 2 5.00 6.00 6.0 7 ▂▃▆▇▅ -0.53 -0.21 2.98
exp_3 0 1.00 5.38 1.26 2 5.00 6.00 6.0 7 ▂▃▆▇▅ -0.63 -0.15 2.69
exp_4 0 1.00 5.31 1.32 1 4.00 6.00 6.0 7 ▁▁▂▃▇ -0.74 0.06 3.25
exp_5 0 1.00 5.10 1.39 1 4.00 5.00 6.0 7 ▁▂▃▅▇ -0.60 -0.12 2.96
exp_6 0 1.00 5.43 1.24 1 5.00 6.00 6.0 7 ▁▁▂▃▇ -0.68 0.02 3.57
int_1 0 1.00 5.32 1.24 1 4.00 6.00 6.0 7 ▁▁▃▃▇ -0.49 -0.35 3.48
int_2 0 1.00 5.33 1.31 1 4.00 6.00 6.0 7 ▁▁▃▃▇ -0.47 -0.50 3.30
int_3 0 1.00 5.23 1.17 1 4.00 5.00 6.0 7 ▁▁▅▅▇ -0.10 -0.78 3.60
int_4 0 1.00 5.32 1.22 1 4.00 6.00 6.0 7 ▁▁▃▃▇ -0.45 -0.37 3.53
ben_1 0 1.00 5.23 1.21 1 4.00 5.00 6.0 7 ▁▁▅▅▇ -0.23 -0.45 3.51
ben_2 0 1.00 5.17 1.26 1 4.00 5.00 6.0 7 ▁▁▆▅▇ -0.22 -0.34 3.31
ben_3 0 1.00 5.27 1.26 1 4.00 5.00 6.0 7 ▁▁▃▃▇ -0.63 0.26 3.37
ben_4 0 1.00 5.00 1.23 1 4.00 5.00 6.0 7 ▁▂▆▅▇ -0.09 -0.56 3.26
meas_rep 0 1.00 1.49 0.50 1 1.00 1.00 2.0 2 ▇▁▁▁▇ 0.05 -2.00 1.02
tsm_1 0 1.00 1.87 0.87 1 1.00 2.00 2.0 4 ▇▇▁▃▁ 0.69 -0.34 2.46
tsm_2 0 1.00 2.05 0.95 1 1.00 2.00 3.0 4 ▇▇▁▅▂ 0.51 -0.73 2.06
tsm_3 0 1.00 2.05 0.97 1 1.00 2.00 3.0 4 ▇▇▁▃▂ 0.62 -0.62 2.01
tsm_4 0 1.00 2.29 1.01 1 1.00 2.00 3.0 4 ▆▇▁▆▃ 0.23 -1.06 1.69
age 21 0.96 22.89 2.95 18 20.00 22.00 25.0 33 ▇▆▅▁▁ 0.72 -0.06 3.43
semester 21 0.96 5.86 3.68 1 3.00 5.00 9.0 19 ▇▅▅▁▁ 0.57 -0.44 3.57
sex 21 0.96 1.71 0.47 1 1.00 2.00 2.0 3 ▃▁▇▁▁ -0.70 -0.91 2.74
tch_1 0 1.00 -107.78 313.71 -999 1.00 2.00 4.0 4 ▁▁▁▁▇ -2.48 4.18 2.84
tch_2 0 1.00 -265.55 443.71 -999 -999.00 1.00 3.5 4 ▃▁▁▁▇ -1.05 -0.91 1.65
tch_3 0 1.00 -278.85 450.45 -999 -999.00 1.00 4.0 4 ▃▁▁▁▇ -0.97 -1.06 1.60
tch_4 0 1.00 -111.51 318.42 -999 1.00 2.00 4.0 4 ▁▁▁▁▇ -2.42 3.89 2.79
tch_5 0 1.00 -276.95 449.52 -999 -999.00 1.00 3.0 4 ▃▁▁▁▇ -0.98 -1.04 1.61
Experitse 0 1.00 5.35 1.12 2 4.67 5.50 6.0 7 ▁▂▅▇▅ -0.59 -0.09 3.00
Integrity 0 1.00 5.30 1.08 2 4.50 5.50 6.0 7 ▁▃▆▇▅ -0.28 -0.60 3.05
Benevolence 0 1.00 5.17 1.06 1 4.25 5.25 6.0 7 ▁▁▆▇▆ -0.16 -0.16 3.93
TSM 0 1.00 2.06 0.67 1 1.50 2.00 2.5 4 ▆▆▇▂▁ 0.33 -0.25 2.88

3.2 METI

3.2.1 Dimensionality (CFA)

First we specified two consecutive threedimensional CFA models

# onedimensional model
cfa_meti_model_1d <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
                             int_1 + int_2 + int_3 + int_4 +
                             ben_1 + ben_2 + ben_3 + ben_4"

cfa1d_meti_1 <- cfa(cfa_meti_model_1d, data = data_rbt%>%filter(meas_rep == 1))
cfa1d_meti_2 <- cfa(cfa_meti_model_1d, data = data_rbt%>%filter(meas_rep == 2))

cfa_meti_model_1 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                     int =~ int_1 + int_2 + int_3 + int_4 
                     ben =~ ben_1 + ben_2 + ben_3 + ben_4
                     int_1 ~~ int_2"

cfa_meti_model_2 <- "exp =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                     int =~ int_1 + int_2 + int_3 + int_4 
                     ben =~ ben_1 + ben_2 + ben_3 + ben_4
                     ben_1 ~~ ben_3
                     ben_1 ~~ ben_2"

cfa_meti_1 <- cfa(cfa_meti_model_1, data = data_rbt%>%filter(meas_rep == 1))

# define a function which prints the fit
fpf <- function(x){  
  fm_tmp <- fitmeasures(x)
  return(cat(sprintf(
          "χ^2^ = %s, _df_ = %s, CFI = %s, TLI = %s, RMSEA = %s, SRMR = %s, SRMR~between~ = %s, SRMR~within~ = %s",
           round(fm_tmp[c("chisq")],3), 
                 fm_tmp[c("df")],
           round(fm_tmp[c("cfi")],3),
           round(fm_tmp[c("tli")],3),
           round(fm_tmp[c("rmsea")],3),
           round(fm_tmp[c("srmr")],3),
           round(fm_tmp[c("srmr_between")],3),
           round(fm_tmp[c("srmr_within")],3))))
}

# print the fit for cfa_meti_1
fpf(cfa1d_meti_1)

χ2 = 588.932, df = 77, CFI = 0.811, TLI = 0.776, RMSEA = 0.157, SRMR = 0.084, SRMRbetween = NA, SRMRwithin = NA

fpf(cfa_meti_1)

χ2 = 194.174, df = 73, CFI = 0.955, TLI = 0.944, RMSEA = 0.078, SRMR = 0.049, SRMRbetween = NA, SRMRwithin = NA

cfa_meti_2 <- cfa(cfa_meti_model_2, data = data_rbt%>%filter(meas_rep == 2))
fpf(cfa1d_meti_2)

χ2 = 936.555, df = 77, CFI = 0.759, TLI = 0.715, RMSEA = 0.208, SRMR = 0.099, SRMRbetween = NA, SRMRwithin = NA

fpf(cfa_meti_2)

χ2 = 212.262, df = 72, CFI = 0.961, TLI = 0.95, RMSEA = 0.087, SRMR = 0.04, SRMRbetween = NA, SRMRwithin = NA

anova(cfa1d_meti_1, cfa_meti_1)
## Chi-Squared Difference Test
## 
##              Df     AIC     BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
## cfa_meti_1   73  9738.4  9853.5 194.17                                  
## cfa1d_meti_1 77 10125.1 10225.9 588.93     394.76       4  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(cfa1d_meti_2, cfa_meti_2)
## Chi-Squared Difference Test
## 
##              Df    AIC    BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
## cfa_meti_2   72 8522.8 8639.9 212.26                                  
## cfa1d_meti_2 77 9237.1 9336.5 936.55     724.29       5  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In an next step we ran a two-level CFA …

# onedimensional model
mcfa1d_meti_model <- "level: 1
                    meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
                            int_1 + int_2 + int_3 + int_4 +
                            ben_1 + ben_2 + ben_3 + ben_4
                    
                    level: 2
                    meti =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6 +
                            int_1 + int_2 + int_3 + int_4 +
                            ben_1 + ben_2 + ben_3 + ben_4"



mcfa_meti_model <- "level: 1
                    exp_w =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                    int_w =~ int_1 + int_2 + int_3 + int_4 
                    ben_w =~ ben_1 + ben_2 + ben_3 + ben_4
                    int_1 ~~ int_2
                    int_3 ~~ ben_4 
                    ben_1 ~~ ben_2 
                 
                    
                    level: 2
                    exp_b =~ exp_1 + exp_2 + exp_3 + exp_4 + exp_5 + exp_6
                    int_b =~ int_1 + int_2 + int_3 + int_4 
                    ben_b =~ ben_1 + ben_2 + ben_3 + ben_4"
mcfa1d_meti <- cfa(mcfa1d_meti_model, data = data_rbt, cluster = "session")
mcfa_meti <- cfa(mcfa_meti_model, data = data_rbt, cluster = "session")
fpf(mcfa1d_meti)

χ2 = 764.777, df = 154, CFI = 0.894, TLI = 0.875, RMSEA = 0.087, SRMR = 0.271, SRMRbetween = 0.146, SRMRwithin = 0.125

fpf(mcfa_meti)

χ2 = 277.759, df = 145, CFI = 0.977, TLI = 0.971, RMSEA = 0.042, SRMR = 0.172, SRMRbetween = 0.091, SRMRwithin = 0.081

anova(mcfa1d_meti, mcfa_meti)
## Chi-Squared Difference Test
## 
##              Df   AIC   BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
## mcfa_meti   145 18147 18484 277.76                                  
## mcfa1d_meti 154 18616 18915 764.78     487.02       9  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.2.2 Reliability (McDonalds \(\omega\))

# First Measurement
## Expertise
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("exp")))$est
## [1] 0.9247211
# Integrity
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("int")))$est
## [1] 0.8806539
# Benevolence
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 1) %>% select(starts_with("ben")))$est
## [1] 0.8513474
# Second Measurement
## Expertise
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("exp")))$est
## [1] 0.9490452
# Integrity
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("int")))$est
## [1] 0.9120384
# Benevolence
MBESS::ci.reliability(data_rbt%>%filter(meas_rep == 2) %>% select(starts_with("ben")))$est
## [1] 0.9003258

3.3 Topic specific multiplism

3.3.1 Dimensionality (CFA)

First we specified two consecutive onedimensional CFA models

cfa_tsm_model <- "tsm =~ a*tsm_1 + a*tsm_2 + a*tsm_3 + a*tsm_4
                  tsm_2 ~~ tsm_4"

cfa_tsm_1 <- cfa(cfa_tsm_model, data = data_rbt%>%filter(meas_rep == 1))
fpf(cfa_tsm_1)

χ2 = 5.302, df = 4, CFI = 0.991, TLI = 0.987, RMSEA = 0.035, SRMR = 0.048, SRMRbetween = NA, SRMRwithin = NA

cfa_tsm_2 <- cfa(cfa_tsm_model, data = data_rbt%>%filter(meas_rep == 2))
fpf(cfa_tsm_2)

χ2 = 6.072, df = 4, CFI = 0.99, TLI = 0.985, RMSEA = 0.045, SRMR = 0.055, SRMRbetween = NA, SRMRwithin = NA

In an next step, we ran a two-level CFA …

mcfa_tsm_model <- "level: 1
                    tsm_w =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
                    tsm_2 ~~ tsm_4
                    

                    level: 2
                    tsm_b =~ tsm_1 + tsm_2 + tsm_3 + tsm_4
                    tsm_2 ~~ 0*tsm_2"

mcfa_tsm <- cfa(mcfa_tsm_model, data = data_rbt, cluster = "session")
fpf(mcfa_tsm)

χ2 = 4.422, df = 4, CFI = 0.999, TLI = 0.996, RMSEA = 0.014, SRMR = 0.109, SRMRbetween = 0.09, SRMRwithin = 0.019

3.3.2 Reliability (McDonalds \(\omega\))

# First Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep ==1) %>% select(starts_with("tsm")))$est
## [1] 0.6509044
# Second Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep ==2) %>% select(starts_with("tsm")))$est
## [1] 0.7125262

3.4 Treatment check

3.4.1 Dimensionality (CFA)

We specified two consecutive onedimensional CFA models

cfa_tch_model <- "tch =~ tch_1 + tch_2 + tch_3 + tch_4 + tch_5
tch_1 ~~ tch_4
tch_3 ~~ tch_5"

cfa_tch_1 <- cfa(cfa_tch_model, 
                 data = data_tch_n%>%
                   filter(meas_rep == 1))
fpf(cfa_tch_1)

χ2 = 2.412, df = 3, CFI = 1, TLI = 1.003, RMSEA = 0, SRMR = 0.006, SRMRbetween = NA, SRMRwithin = NA

cfa_tch_2 <- cfa(cfa_tch_model, data = data_tch_n%>%filter(meas_rep == 2))
fpf(cfa_tch_2)

χ2 = 6.538, df = 3, CFI = 0.997, TLI = 0.991, RMSEA = 0.083, SRMR = 0.005, SRMRbetween = NA, SRMRwithin = NA

3.4.2 Reliability (Ordinal McDonalds \(\omega\))

# First Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep == 2) %>% select(starts_with("tch")))$est
## [1] 0.8625644
# Second Measurement
MBESS::ci.reliability(data_rbt %>% filter(meas_rep == 2) %>% select(starts_with("tch")))$est
## [1] 0.8625644

3.5 Table of (M)CFA fit-indices

tibble(`1d CFA METI 1` = fitmeasures(cfa1d_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA METI 2` = fitmeasures(cfa1d_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `3d CFA METI 1` = fitmeasures(cfa_meti_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `3d CFA METI 2` = fitmeasures(cfa_meti_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d MCFA METI` = fitmeasures(mcfa1d_meti)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `3d MCFA METI` = fitmeasures(mcfa_meti)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TSM 1` = fitmeasures(cfa_tsm_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TSM 2` = fitmeasures(cfa_tsm_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d MCFA TSM` = fitmeasures(mcfa_tsm)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TCH 1` = fitmeasures(cfa_tch_1)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       `1d CFA TCH 2` = fitmeasures(cfa_tch_2)[c("chisq", "df", "cfi", "tli", "rmsea", 
                                              "srmr", "srmr_between", "srmr_within", "bic", "aic")],
       rownames = c("χ^2^", "_df_", "CFI", "TLI", "RMSEA", "SRMR", "SRMR~between~", "SRMR~within~", "BIC", "AIC"))%>%
  column_to_rownames(var = "rownames")%>%
  knitr::kable(., digits = 3)
1d CFA METI 1 1d CFA METI 2 3d CFA METI 1 3d CFA METI 2 1d MCFA METI 3d MCFA METI 1d CFA TSM 1 1d CFA TSM 2 1d MCFA TSM 1d CFA TCH 1 1d CFA TCH 2
χ2 588.932 936.555 194.174 212.262 764.777 277.759 5.302 6.072 4.422 2.412 6.538
df 77.000 77.000 73.000 72.000 154.000 145.000 4.000 4.000 4.000 3.000 3.000
CFI 0.811 0.759 0.955 0.961 0.894 0.977 0.991 0.990 0.999 1.000 0.997
TLI 0.776 0.715 0.944 0.950 0.875 0.971 0.987 0.985 0.996 1.003 0.991
RMSEA 0.157 0.208 0.078 0.087 0.087 0.042 0.035 0.045 0.014 0.000 0.083
SRMR 0.084 0.099 0.049 0.040 0.271 0.172 0.048 0.055 0.109 0.006 0.005
SRMRbetween NA NA NA NA 0.146 0.091 NA NA 0.090 NA NA
SRMRwithin NA NA NA NA 0.125 0.081 NA NA 0.019 NA NA
BIC 10225.876 9336.494 9853.512 8639.946 18914.759 18484.147 2791.127 2653.377 5445.587 1547.731 1750.473
AIC 10125.120 9237.119 9738.362 8522.826 18616.055 18147.038 2769.537 2632.082 5360.243 1512.868 1712.703

4 Results of the treatmentcheck

4.1 Plot

res_tch_data <- data_tch%>%
  gather(variable, value, starts_with("tch_"))%>%
  group_by(treat, variable, value)%>%
  mutate(value = as.character(value)) %>% 
  summarize(freq = n())%>%
  ungroup()%>%
  mutate(treat = case_when(treat == "GB" ~ "Greyed out badges",
                           treat == "CB" ~ "Colored badges",
                           treat == "CC" ~ "Control condition",
                           T ~ treat),
         value = case_when(value =="-999" ~ "don't know",
                           T ~ value),
         variable = case_when(variable == "tch_1" ~ "item 1", 
                              variable == "tch_2" ~ "item 2", 
                              variable == "tch_3" ~ "item 3", 
                              variable == "tch_4" ~ "item 4", 
                              variable == "tch_5" ~ "item 5"),
         Frequency = freq)
## `summarise()` regrouping output by 'treat', 'variable' (override with `.groups` argument)
res_tch_plot <- ggplot(res_tch_data, aes(variable, value)) + 
  geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
  scale_size_continuous(range = c(3,15)) + 
  scale_color_gradient(low = "grey95", high = "grey65") +
  guides(color=guide_legend(), size = guide_legend()) +
  facet_wrap(~treat) + 
  theme_ipsum_ps() + 
  ggtitle("Results of the treatment check", "Frequency per item and experimental condition") + 
  ylab("") + 
  xlab("")
  
res_tch_plot

res_tch_plot_pub <- ggplot(res_tch_data %>%
                         mutate(treat = case_when(treat == "Greyed out badges" ~ "Grey badges",
                                                  TRUE ~ treat)), 
                                aes(variable, value)) + 
 geom_point(aes(size = Frequency, color = Frequency), shape = 15) +
 scale_size_continuous(range = c(3,15)) + 
 scale_color_gradient(low = "grey95", high = "grey65") +
 guides(color=guide_legend(), size = guide_legend()) +
 facet_wrap(~treat) + 
 theme_ipsum_ps() + 
 ylab("") + 
 xlab("")


#ggsave("8_reproducible_documentation_of_analyses/Fig/Fig2.jpg", res_tch_plot_pub, width = 130, height = 50, units = "mm", dpi = 300, scale = 2.4)

4.2 Effect size

res_tch_data_A <- data_tch_n%>%
  filter(treat != "CC")%>%
  filter(tch_1 != -999)

effsize::VD.A(tch_1 ~ treat, data = res_tch_data_A)   
## 
## Vargha and Delaney A
## 
## A estimate: 0.837694 (large)

5 Graphical exploration

5.1 Plot Hyp 1

data_scales_lables%>%
  gather(Variable, Value, Experitse, Integrity, Benevolence)%>%
  ggplot(., aes(x = Treatment, y = Value)) + 
  geom_violin(adjust = 1.5) +
  stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
  coord_flip() +
  facet_wrap(~ Variable, nrow = 1) + 
  labs(title = "Graphical overview (Hyp 1)",
       subtitle = "Violinplots and means ± 1*SD",
       caption = "") +
  ylim(1,7) +
  hrbrthemes::theme_ipsum_ps()

# Export data for combined plot
#data_scales_lables %>% 
#  mutate(Study = "Study 1") %>% 
#  write_csv(., "8_reproducible_documentation_of_analyses/Fig/data_scales_lables_study_1.csv")

5.2 Descriptive Effect Sizes Hyp 1

A_GB_CC <- data_scales_lables%>%
  filter(treat != "CB")%>%
  mutate(treat = as.character(treat))%>%
  effsize::VD.A(Integrity ~ treat, data = .)%>%
  .$estimate


A_CC_CB <- data_scales_lables%>%
  filter(treat != "GB")%>%
  mutate(treat = as.character(treat))%>%
  effsize::VD.A(Integrity ~ treat, data = .)%>%
  .$estimate


A_GB_CB <- data_scales_lables%>%
  filter(treat != "CC")%>%
  mutate(treat = as.character(treat))%>%
  effsize::VD.A(Integrity ~ treat, data = .)%>%
  .$estimate
GB CC
CC A = 0.59, d = 0.32
CB A = 0.66, d = 0.57 A = 0.58, d = 0.29

5.3 Hyp 2/3

plot_hyp2_1 <- ggplot(data_scales_lables, 
                      aes(x=`Topic specific multiplism`, y = Integrity)) + 
  geom_jitter() +
  facet_wrap(~ Treatment, nrow = 1) + 
  labs(title = "Graphical overview (Hyp 2/3)",
       subtitle = "Jitter plot per treatment") +
  hrbrthemes::theme_ipsum_ps()


plot_hyp2_1 + stat_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

plot_hyp2_1 + stat_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'

5.4 Descriptive Effect Sizes Hyp 3/4

Spearman and Kendall correlations:

r_GB <- round(cor(data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(`Topic specific multiplism`),
                      method = "pearson", use = "pairwise.complete.obs"), 2)
t_GB <- round(cor(data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "GB")%>%
                        pull(`Topic specific multiplism`),
                      method = "kendall", use = "pairwise.complete.obs"), 2)

r_CC <- round(cor(data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(`Topic specific multiplism`),
                      method = "pearson", use = "pairwise.complete.obs"), 2)
t_CC <- round(cor(data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CC")%>%
                        pull(`Topic specific multiplism`),
                      method = "kendall", use = "pairwise.complete.obs"), 2)

r_CB <- round(cor(data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(`Topic specific multiplism`),
                      method = "pearson", use = "pairwise.complete.obs"), 2)
t_CB <- round(cor(data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(Integrity),
                      data_scales_lables%>%
                        filter(treat == "CB")%>%
                        pull(`Topic specific multiplism`),
                      method = "kendall", use = "pairwise.complete.obs"), 2)
GB CC CB
\(r(integrity, multiplism)\) -0.46 -0.36 -0.35
\(\tau(integrity, multiplism)\) -0.33 -0.24 -0.26

5.5 Hyp 4

data_scales_lables%>%
  mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
                               treat == "CB" ~ "Colored badges",
                               treat == "CC" ~ "Control condition",
                               T ~ "treat"))%>%
  ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) + 
  geom_violin(adjust = 1.5, alpha = .5) +
  stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
  coord_flip() +
  labs(title = "Graphical overview (Hyp 4)",
       subtitle = "Violinplots and means ± 1*SD",
       caption = "") +
  xlab("") +
  ylim(1,4) +
  hrbrthemes::theme_ipsum_ps()

 fig4 <- data_scales_lables%>%
   mutate(Treatment = case_when(treat == "GB" ~ "Greyed out badges",
                                treat == "CB" ~ "Colored badges",
                                treat == "CC" ~ "Control condition",
                                T ~ "treat"))%>%
   ggplot(., aes(x = Treatment, y = `Topic specific multiplism`)) + 
   geom_violin(adjust = 1.5, alpha = .5) +
   stat_summary(fun.data = "mean_sdl", fun.args = list(mult = 1)) +
   coord_flip() +
   xlab("") +
   ylim(1,4) +
   hrbrthemes::theme_ipsum_ps()

5.6 Descriptive Effect Sizes Hyp 4

A_mult_GB_CC <- effsize::VD.A(`Topic specific multiplism` ~ treat, 
                               data = data_scales_lables%>%
                                 filter(treat != "CB")%>%
                                 mutate(treat = as.character(treat)))$estimate
d_mult_GB_CC <- qnorm(A_mult_GB_CC)*sqrt(2)


A_mult_CC_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat, 
                               data = data_scales_lables%>%
                                 filter(treat != "GB")%>%
                                 mutate(treat = as.character(treat)))$estimate
d_mult_CC_CB <- qnorm(A_mult_CC_CB)*sqrt(2)


A_mult_GB_CB <- effsize::VD.A(`Topic specific multiplism` ~ treat, 
                               data = data_scales_lables%>%
                                 filter(treat != "CC")%>%
                                 mutate(treat = as.character(treat)))$estimate
d_mult_GB_CB <- qnorm(A_mult_GB_CB)*sqrt(2)
GB CC
CC A = 0.43, d = -0.26
CB A = 0.43, d = -0.25 A = 0.5, d = 0.01

6 Inference statistics (Hyp 1)

6.1 Description of the variables

All of the following analyses are based on the data frame object data_scales_wide which is why we describe it here somewhat more detailed.
All analyses are based on measured variables Integrity (Integrity, is information source sincere, honest, just, unselfish and fair?) and Topic Specific Multiplism (TSM, is knowledge and knowing about a topic arbitrary?). As this data set is in wide format, the experimental conditions are encoded wtihin the variable names:

  • GB means, that the participants of the study could learn from the grey out badges (ans corresponding explanations) within the abstract, that the authors of the study denied to use Open Practices
  • CC means, that the participants of the study could not learn if or if not the authors of the study used Open Practices
  • CBmeans, that the participants of the study could learn from the colored badges (ans corresponding explanations) within the abstract, that the authors of the study used Open Practices

Finally, the variable session identified the study participants.

If we look descriptively at these variables:

data_scales_wide%>%
  my_skim(.)
Data summary
Name Piped data
Number of rows 270
Number of columns 13
_______________________
Column type frequency:
character 1
numeric 12
________________________
Group variables None

Variable type: character

skim_variable n_missing complete_rate min max empty n_unique whitespace
session 0 1 64 64 0 270 0

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist skewn kurto maxabszvalue
Benevolence_CB 100 0.63 5.44 0.99 2.75 4.75 5.25 6.25 7 ▁▂▇▅▆ -0.08 -0.71 2.71
Benevolence_CC 97 0.64 5.20 0.96 2.75 4.50 5.25 6.00 7 ▁▅▇▆▃ 0.03 -0.61 2.55
Benevolence_GB 86 0.68 4.89 1.15 1.00 4.00 4.75 5.75 7 ▁▁▇▆▅ -0.14 -0.02 3.39
Experitse_CB 100 0.63 5.63 1.10 2.33 4.88 5.83 6.50 7 ▁▂▃▇▇ -0.78 0.13 2.99
Experitse_CC 97 0.64 5.35 0.97 2.50 4.67 5.33 6.00 7 ▁▃▆▇▅ -0.25 -0.48 2.94
Experitse_GB 86 0.68 5.09 1.19 2.00 4.29 5.33 6.00 7 ▂▂▅▇▃ -0.59 -0.32 2.59
Integrity_CB 100 0.63 5.62 0.98 3.25 5.00 5.75 6.50 7 ▂▅▇▇▇ -0.31 -0.83 2.42
Integrity_CC 97 0.64 5.34 0.97 2.75 4.50 5.50 6.00 7 ▁▃▇▇▅ -0.22 -0.62 2.66
Integrity_GB 86 0.68 4.97 1.18 2.00 4.00 5.00 6.00 7 ▂▅▆▇▃ -0.09 -0.76 2.53
TSM_CB 100 0.63 2.01 0.67 1.00 1.50 2.00 2.50 4 ▇▇▇▂▁ 0.38 -0.34 2.95
TSM_CC 97 0.64 1.99 0.61 1.00 1.50 2.00 2.50 4 ▆▆▇▁▁ 0.25 -0.12 3.30
TSM_GB 86 0.68 2.18 0.72 1.00 1.75 2.25 2.56 4 ▅▆▇▂▂ 0.24 -0.46 2.54

6.2 Investigating the missingness

6.2.1 Missingness per Variable

library(naniar)
## 
## Attaching package: 'naniar'
## The following object is masked from 'package:skimr':
## 
##     n_complete
visdat::vis_miss(data_scales_wide%>%select(-session)) + 
  ggtitle("Missingness per Variable") + 
  theme(plot.margin = margin(1, 2, 1, 1, "cm"))

6.2.2 Marginal distributions Integrity_GB

library(patchwork)
ggplot(data_scales_wide, aes(Integrity_GB, Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       Integrity_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       Integrity_GB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.3 Marginal distributions Integrity_CC

ggplot(data_scales_wide, aes(Integrity_GB, Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       Integrity_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       Integrity_CC)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.4 Marginal distributions Integrity_CB

ggplot(data_scales_wide, aes(Integrity_GB, Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       Integrity_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       Integrity_CB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.5 Marginal distributions TSM_GB

ggplot(data_scales_wide, aes(Integrity_GB, TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       TSM_GB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       TSM_GB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.6 Marginal distributions TSM_CC

ggplot(data_scales_wide, aes(Integrity_GB, TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       TSM_CC)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       TSM_CC)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.7 Marginal distributions TSM_CB

ggplot(data_scales_wide, aes(Integrity_GB, TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CC, TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(Integrity_CB, TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_GB,       TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CC,       TSM_CB)) + 
  geom_miss_point(alpha = .2) +
  theme_ipsum_ps() +
ggplot(data_scales_wide, aes(TSM_CB,       TSM_CB)) + 
  geom_miss_point(alpha = .2) + plot_layout(guides = 'collect') +
  theme_ipsum_ps()

6.2.8 Cohen’s d of missing/not-missing per variable

d_matrix_missings <- miceadds::mi_dstat(data_scales_wide%>%select(starts_with("Int"), starts_with("TSM")))%>%
  round(.,4)

knitr::kable(d_matrix_missings, caption = "Boxplot of Cohen's d of missing/not-missing per variable")
Boxplot of Cohen’s d of missing/not-missing per variable
Integrity_CB Integrity_CC Integrity_GB TSM_CB TSM_CC TSM_GB
Integrity_CB NaN -0.2103 -0.0096 NaN -0.1009 -0.3312
Integrity_CC -0.1596 NaN -0.0074 0.1056 NaN 0.3622
Integrity_GB 0.1045 0.1696 NaN -0.1339 0.1018 NaN
TSM_CB NaN -0.2103 -0.0096 NaN -0.1009 -0.3312
TSM_CC -0.1596 NaN -0.0074 0.1056 NaN 0.3622
TSM_GB 0.1045 0.1696 NaN -0.1339 0.1018 NaN 
boxplot(as.vector(d_matrix_missings))
title("Boxplot of Cohen's d of missing/not-missing per variable")

6.3 Imputation

M <- 1000
out <- mice(data = data_scales_wide%>%select(-session),
            m = M,
            meth=c("norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm"), 
            diagnostics = FALSE,
            printFlag = F,
            seed = 83851)

6.4 Check of first 10 imputation

out_first10 <-   mice(data = data_scales_wide%>%select(-session),
            m = 10,
            meth=c("norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm",
                   "norm","norm","norm"), 
            diagnostics = FALSE,
            printFlag = F,
            seed = 83851)

densityplot(out_first10)

plot(out_first10)

6.5 Parameter and BF estimation

# Set up the matrices for the estimates ##############
# setup of matrices to store multiple estimates
mulest_hyp1 <- matrix(0,nrow=M,ncol=3) 
# and covariance matrices
covwithin_hyp1 <- matrix(0,nrow=3,ncol=3) 


# Estimate the coefficients for each data frame ######
for(i in 1:M) {
 within_hyp1 <- lm(cbind(Integrity_GB,Integrity_CC,Integrity_CB)~1, 
                   data=mice::complete(out,i)) # estimate the means of the three variables
 mulest_hyp1[i,]<-coef(within_hyp1)[1:3] # store these means in the matrix `mulres`
 covwithin_hyp1<-covwithin_hyp1 + 1/M * vcov(within_hyp1)[1:3,1:3] # compute the averages 
}


# Compute the average of the estimates ###############
estimates_hyp1 <- colMeans(mulest_hyp1)
names(estimates_hyp1) <- c("Integrity_GB","Integrity_CC","Integrity_CB")
covbetween_hyp1 <- cov(mulest_hyp1) # between covariance matrix
covariance_hyp1 <- covwithin_hyp1 + (1+1/M)*covbetween_hyp1 # total variance


# Determine the effective and real sample sizes ######
samp_hyp1 <- nrow(data_scales_wide) # real sample size
nucom_hyp1<-samp_hyp1-length(estimates_hyp1)

# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp1 <- (1+1/M)*(1/length(estimates_hyp1))* 
  sum(diag(covbetween_hyp1 %*% 
             MASS::ginv(covariance_hyp1))) # ... (43)
nuold_hyp1<-(M-1)/(lam_hyp1^2) # ... (44)
nuobs_hyp1<-(nucom_hyp1+1)/(nucom_hyp1+3)*nucom_hyp1*(1-lam_hyp1) # ... (46)
nu_hyp1<- nuold_hyp1*nuobs_hyp1/(nuold_hyp1+nuobs_hyp1) # ... (47)
fracmis_hyp1 <- (nu_hyp1+1)/(nu_hyp1+3)*lam_hyp1 + 2/(nu_hyp1+3) # ... (48)
neff_hyp1<-samp_hyp1-samp_hyp1*fracmis_hyp1 # = 172 approx. 2/3* 270

# coerce `covariance` to a list
covariance_hyp1<-list(covariance_hyp1)


# Test the hypotheses with bain ######################
set.seed(6346)
results_hyp1 <- bain(estimates_hyp1,
               "Integrity_GB<Integrity_CC<Integrity_CB;
                Integrity_GB=Integrity_CC=Integrity_CB;
                Integrity_GB<Integrity_CC=Integrity_CB",
                n = neff_hyp1, Sigma=covariance_hyp1,    
                group_parameters=3, joint_parameters = 0)
print(results_hyp1)
## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit   Com   BF.u  BF.c     PMPa  PMPb 
## H1 0.999 0.177 5.656 3878.238 0.970 0.828
## H2 0.000 0.279 0.000 0.000    0.000 0.000
## H3 0.045 0.257 0.174 0.174    0.030 0.025
## Hu                                  0.146
## 
## Hypotheses:
##   H1: Integrity_GB<Integrity_CC<Integrity_CB
##   H2: Integrity_GB=Integrity_CC=Integrity_CB
##   H3: Integrity_GB<Integrity_CC=Integrity_CB
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
summary(results_hyp1)
##      Parameter        n Estimate       lb       ub
## 1 Integrity_GB 168.9123 4.992512 4.830189 5.154835
## 2 Integrity_CC 168.9123 5.327959 5.189725 5.466193
## 3 Integrity_CB 168.9123 5.585524 5.443602 5.727446
results_hyp1$BFmatrix
##              H1       H2           H3
## H1 1.000000e+00 61491369 3.254405e+01
## H2 1.626244e-08        1 5.292458e-07
## H3 3.072758e-02  1889481 1.000000e+00

7 Inference statistics (Hyp 2)

7.1 Parameter estimation

path_mod <- "Integrity_GB ~ TSM_GB
             Integrity_CC ~ TSM_CC             
             Integrity_CB ~ TSM_CB"

# Set up the matrices for the estimates ##############
best_hyp2 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp2 <- matrix(0,nrow=3,ncol=3) # and covariance matrices


# Estimate the coefficients for each data frame ######
for(i in 1:M) {
 path_fit <- sem(path_mod, 
                 data = mice::complete(out, i), 
                 std.lv = T) # estimate the path coefficients 
 best_hyp2[i,] <- parameterestimates(path_fit, standardized = T)%>% # store path coefficients
                     filter(op == "~")%>%
                     pull(std.all) 
 covwithin_hyp2 <- covwithin_hyp2 + # compute the average of the covariance matrices
                  1/M * lavInspect(path_fit, 
                                   "vcov.std.all")[c("Integrity_GB~TSM_GB", 
                                                     "Integrity_CC~TSM_CC", 
                                                     "Integrity_CB~TSM_CB"),
                                                   c("Integrity_GB~TSM_GB", 
                                                     "Integrity_CC~TSM_CC", 
                                                     "Integrity_CB~TSM_CB")]
}

# Compute the average of the estimates ###############
estimates_hyp2 <- colMeans(best_hyp2)
names(estimates_hyp2) <- c("Int_on_TSM_GB", "Int_on_TSM_CC", "Int_on_TSM_CB")
round(estimates_hyp2, 2)
## Int_on_TSM_GB Int_on_TSM_CC Int_on_TSM_CB 
##         -0.44         -0.41         -0.32

7.2 Visual path model

\usetikzlibrary{arrows,positioning}
\usetikzlibrary{calc}
\begin{tikzpicture}[auto,>=latex,align=center,
latent/.style={circle,draw, thick,inner sep=0pt,minimum size=10mm},
manifest/.style={rectangle,draw, thick,inner sep=2pt,minimum size=15mm},
manifestcov/.style={rectangle,draw, thick,inner sep=1pt,minimum size=8mm},
mean/.style={regular polygon,regular polygon sides=3,draw,thick,inner sep=0pt,minimum size=10mm},
paths/.style={->, thick, >=latex, font = \footnotesize\sffamily, fill = white},
variance/.style={<->, thin, >=latex, bend left=270, looseness=2.5, font = \footnotesize},
covariance/.style={<->, thin, >=latex, bend left=290, looseness=.5, font = \footnotesize},
dotsdots/.style={circle, draw, fill = black, minimum size = 1mm, maximum size = 1mm}
]

%% AVs
\node [manifest] (AV1) at (0, -2)       {$int_{gb}$};
\node [manifest] (AV2) [right = of AV1] {$int_{cc}$};
\node [manifest] (AV3) [right = of AV2] {$int_{cb}$};


%% UVs
\node [manifest] (UV1) at (0, 2)        {$tsm_{gb}$};    
\node [manifest] (UV2) [right =of UV1]  {$tsm_{cc}$};    
\node [manifest] (UV3) [right =of UV2]  {$tsm_{cb}$};    

%% Paths
\draw [paths] (UV1) to node [midway]{-.44} (AV1);
\draw [paths] (UV2) to node [midway]{-.41} (AV2);
\draw [paths] (UV3) to node [midway]{-.32} (AV3);

%\draw [covariance] (UV2) to node {} (UV1);
\end{tikzpicture}

covbetween_hyp2 <- cov(best_hyp2) # between covariance matrix
covariance_hyp2 <- covwithin_hyp2 + (1+1/M)*covbetween_hyp2 # total variance

# Determine the effective and real sample sizes ######
samp_hyp2 <- nrow(data_scales_wide) # real sample size
nucom_hyp2 <- samp_hyp2-length(estimates_hyp2)

# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp2 <- (1+1/M)*(1/length(estimates_hyp2))* 
  sum(diag(covbetween_hyp2 %*% 
             MASS::ginv(covariance_hyp2))) # ... (43)
nuold_hyp2 <- (M-1)/(lam_hyp2^2) # ... (44)
nuobs_hyp2 <- (nucom_hyp2+1)/(nucom_hyp2+3)*nucom_hyp2*(1-lam_hyp2) # ... (46)
nu_hyp2 <- nuold_hyp2*nuobs_hyp2/(nuold_hyp2+nuobs_hyp2) # ... (47)
fracmis_hyp2 <- (nu_hyp2+1)/(nu_hyp2+3)*lam_hyp2 + 2/(nu_hyp2+3) # ... (48)
neff_hyp2 <- samp_hyp2-samp_hyp2*fracmis_hyp2 # = 114 approx. 2/3* 270

# coerce `covariance` to a list
covariance_hyp2 <- list(covariance_hyp2)

7.3 Bayes Factor estimation (Hyp 2)

set.seed(6346)
results_hyp2 <- bain(estimates_hyp2, 
                     "Int_on_TSM_GB < 0 & Int_on_TSM_CC < 0 & Int_on_TSM_CB < 0;
                     Int_on_TSM_GB = 0 & Int_on_TSM_CC = 0 & Int_on_TSM_CB = 0",
                    Sigma = covariance_hyp2,
                    n = neff_hyp2,
                    group_parameters = 3,
                    joint_parameters = 0,
                    standardize = F)
print(results_hyp2)
## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit   Com   BF.u  BF.c         PMPa  PMPb 
## H1 1.000 0.161 6.225 10465021.465 1.000 0.862
## H2 0.000 1.281 0.000 0.000        0.000 0.000
## Hu                                      0.138
## 
## Hypotheses:
##   H1: Int_on_TSM_GB<0&Int_on_TSM_CC<0&Int_on_TSM_CB<0
##   H2: Int_on_TSM_GB=0&Int_on_TSM_CC=0&Int_on_TSM_CB=0
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp2$BFmatrix
##              H1           H2
## H1 1.000000e+00 1.273079e+21
## H2 7.854974e-22 1.000000e+00

8 Inference statistics (Hyp 3)

set.seed(6346)
results_hyp3 <- bain(estimates_hyp2, 
                    "Int_on_TSM_GB < Int_on_TSM_CC < Int_on_TSM_CB;
                     (Int_on_TSM_GB, Int_on_TSM_CC) < Int_on_TSM_CB;
                     Int_on_TSM_GB = Int_on_TSM_CC = Int_on_TSM_CB",
                    Sigma = covariance_hyp2,
                    n = neff_hyp2,
                    group_parameters = 3,
                    joint_parameters = 0,
                    standardize = T)

print(results_hyp3)
## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit    Com   BF.u   BF.c   PMPa  PMPb 
## H1 0.547  0.174 3.142  5.733  0.127 0.122
## H2 0.815  0.343 2.375  8.431  0.096 0.092
## H3 10.003 0.518 19.310 19.310 0.778 0.748
## Hu                                  0.039
## 
## Hypotheses:
##   H1: Int_on_TSM_GB<Int_on_TSM_CC<Int_on_TSM_CB
##   H2: (Int_on_TSM_GB,Int_on_TSM_CC)<Int_on_TSM_CB
##   H3: Int_on_TSM_GB=Int_on_TSM_CC=Int_on_TSM_CB
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp3$BFmatrix
##           H1       H2        H3
## H1 1.0000000 1.322887 0.1627242
## H2 0.7559224 1.000000 0.1230069
## H3 6.1453658 8.129625 1.0000000

9 Inference statistics (Hyp 4)

# Set up the matrices for the estimates ##############
mulest_hyp4 <- matrix(0,nrow=M,ncol=3) # setup of matrices to store multiple estimates
covwithin_hyp4 <- matrix(0,nrow=3,ncol=3) # and covariance matrices


# Estimate the coefficients for each data frame ######
for(i in 1:M) {
 within_hyp4 <- lm(cbind( TSM_GB, TSM_CC, TSM_CB) ~ 1, 
                   data=mice::complete(out,i)) # estimate the means of the three variables
 mulest_hyp4[i,]<-coef(within_hyp4)[1:3] # store these means in the matrix `mulres`
 covwithin_hyp4<-covwithin_hyp4 + 1/M * vcov(within_hyp4)[1:3,1:3] # compute the averages 
}


# Compute the average of the estimates ###############
estimates_hyp4 <- colMeans(mulest_hyp4)
names(estimates_hyp4) <- c("TSM_GB","TSM_CC","TSM_CB")
covbetween_hyp4 <- cov(mulest_hyp4) # between covariance matrix
covariance_hyp4 <- covwithin_hyp4 + (1+1/M)*covbetween_hyp4 # total variance


# Determine the effective and real sample sizes ######
samp_hyp4 <- nrow(data_scales_wide) # real sample size
nucom_hyp4<-samp_hyp4-length(estimates_hyp4)

# corresponds to Equation (X) in Hoijtink, Gu, Mulder, & Rosseel (2019) ...
lam_hyp4 <- (1+1/M)*(1/length(estimates_hyp4))* 
  sum(diag(covbetween_hyp4 %*% 
             MASS::ginv(covariance_hyp4))) # ... (43)
nuold_hyp4<-(M-1)/(lam_hyp4^2) # ... (44)
nuobs_hyp4<-(nucom_hyp4+1)/(nucom_hyp4+3)*nucom_hyp4*(1-lam_hyp4) # ... (46)
nu_hyp4<- nuold_hyp4*nuobs_hyp4/(nuold_hyp4+nuobs_hyp4) # ... (47)
fracmis_hyp4 <- (nu_hyp4+1)/(nu_hyp4+3)*lam_hyp4 + 2/(nu_hyp4+3) # ... (48)
neff_hyp4<-samp_hyp4-samp_hyp4*fracmis_hyp4 # = 172 approx. 2/3* 270

# coerce `covariance` to a list
covariance_hyp4<-list(covariance_hyp4)


# Test the hypotheses with bain ######################
set.seed(6346)
results_hyp4 <- bain(estimates_hyp4,
               "TSM_GB>TSM_CC>TSM_CB;
                TSM_GB=TSM_CC=TSM_CB;
                (TSM_GB,TSM_CC)>TSM_CB",
                n = neff_hyp4, Sigma=covariance_hyp4,    
                group_parameters=3, joint_parameters = 0)

print(results_hyp4)
## Bayesian informative hypothesis testing for an object of class numeric:
## 
##    Fit   Com   BF.u  BF.c  PMPa  PMPb 
## H1 0.650 0.180 3.602 8.433 0.589 0.507
## H2 0.297 0.510 0.583 0.583 0.095 0.082
## H3 0.658 0.342 1.926 3.708 0.315 0.271
## Hu                               0.141
## 
## Hypotheses:
##   H1: TSM_GB>TSM_CC>TSM_CB
##   H2: TSM_GB=TSM_CC=TSM_CB
##   H3: (TSM_GB,TSM_CC)>TSM_CB
## 
## Note: BF.u denotes the Bayes factor of the hypothesis at hand versus the unconstrained hypothesis Hu. BF.c denotes the Bayes factor of the hypothesis at hand versus its complement.
results_hyp4$BFmatrix
##           H1       H2      H3
## H1 1.0000000 6.181588 1.87024
## H2 0.1617707 1.000000 0.30255
## H3 0.5346909 3.305238 1.00000